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The unsteady flow over a hump model with zero efflux oscillatory flow control is modeled computationally 
using the unsteady Reynolds-averaged Navier-Stokes equations. Three different turbulence models produce 
similar results, and do a reasonably good job predicting the general character of the unsteady surface pressure 
coefficients during the forced cycle. However, the turbulent shear stresses are underpredicted in magnitude 
inside the separation bubble, and the computed results predict too large a (mean) separation bubble compared 
with experiment. These missed predictions are consistent with earlier steady-state results using no-flow-control 
and steady suction, from a 2004 CFD validation workshop for synthetic jets. 


I. Introduction 

Synthetic (or zero efflux) jets are known to be an effective method for controlling separation. * 1 It is believed 
that large unsteady spanwise structures induce momentum transfer across the boundary layer, which in turn reduces 
separation in the mean. As a practical tool, synthetic jets are generally more attractive than control systems using 
steady blowing or suction because complex ducting systems are not required for synthetic jets. Often, unsteady jet 
actuation can be achieved via relatively small piezo-electric plates or electro-mechanical piston arrangements in a 
cavity or cavities located just beneath the body surface. 

Computational methods will play an important role in the design of effective synthetic jet systems. However, 
separated turbulent flows are often difficult to predict with Reynolds-averaged Navier-Stokes (RANS) CFD codes, and 
the complex physics of unsteady flow only serves to make the problems even more difficult to model. It is therefore 
important to perform detailed analyses of the current capabilities of state-of-the-art CFD methods applied to synthetic 
jet problems, and to attempt identify specific deficiencies in the modeling. 

To this end, well-documented experiments on relatively simple configurations (or “unit problems”) are particularly 
useful. A CFD validation workshop for synthetic jets and turbulent separation control (CFDVAL2004) was held in 
2004, and included three such unit problems. 2,3 One of the cases was the flow over a two-dimensional wall-mounted 
hump-shaped body, termed the ”hump model” (it was actually mounted on a splitter plate). During the CFDVAL2004 
workshop, the effort was focused primarily on steady-state hump model results using no-flow-control and steady 
suction flow control. There was only limited comparison with oscillatory synthetic jet data because the unsteady 
experiment was not completed at the time of the workshop. The steady experimental results were documented in 
Greenblatt et al., 4 with wall shear stress measurements for the no-flow-control case documented in Naughton et al. 5 
After the CFDVAL2004 workshop, extensive time-accurate data became available for the oscillatory flow-control 
case. 6 Measurement uncertainties were also documented (±0.001 for surface pressure coefficients, less than 3% for 
velocities, and ±10% for turbulence quantities). Overall, the hump model case has generated a great deal of interest 
and has spawned many CFD studies. 7-14 

Rumsey et al. 2 provided a summary of the hump model CFDVAL2004 workshop results compared with experi- 
ment. There were two key points made concerning this case. First, the side plates used in the experiment to improve 
two-dimensionality were found to cause blockage that, if not modeled in CFD, resulted in noticeable overprediction 
of pressures over most of the hump surface. Second, RANS CFD results consistently predicted too long a separation 
bubble for both the no-flow-control and the steady-suction cases. This overprediction was surmised to be due to the 
fact that the modeled turbulent shear stresses in the separated region were too low in magnitude (by as much as a 
factor of 2 - 3), compared with experimental results. More computationally intensive methods that involved large 
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eddy simulation (LES) (such as blended RANS-LES techniques) or direct numerical simulation (DNS) appeared to 
hold promise for predicting improved separation lengths, but results were inconclusive at the time of the workshop. 

The current study is an application of unsteady Reynolds-averaged Navier-Stokes (URANS) to the oscillatory case 
using the CFL3D code 15 with three different turbulence models: the Spalart-Allmaras (SA) model, 16 Menter’s k-w 
SST model, 17 and the explicit algebraic stress model in k-u form (EASM-ko). 18 The jet slot and cavity are included 
in the computations. In a related effort, Balakumar 12 performed URANS computations with higher order spatial 
accuracy using the SST model, and found that the computations did not predict the correct mean separation length 
(predictions were more than 20% too long). Capizzano et al. 14 also used URANS with the SST model, and their 
mean separation length was similarly overpredicted, by roughly 15%. Computed turbulent shear stress levels were 
significantly underpredicted in magnitude in the separated region. Both Balakumar and Capizzano et al.employed an 
unsteady wall-surface boundary condition rather than including the jet slot and cavity in their computations. Morgan 
et al. 13 used LES (computed at a lower Reynolds number than experiment), including the jet slot and cavity, to achieve 
turbulent shear stress levels in the separated region in generally good agreement with the synthetic jet experiment. 
Mean predicted reattachment length was improved, but was still somewhat longer than experiment (by 10%). 

II. Solution Methodology 

The computer code CFL3D 15 solves the three-dimensional, time-dependent, Reynolds-averaged compressible 
Navier-Stokes equations with an upwind finite-volume formulation (it can also be exercised in two-dimensional mode 
of operation for 2-D cases). Because the code is compressible, the Navier-Stokes equations are averaged using Favre 
averaged variables; e.g., pui/p. CFL3D can solve flows over multiple-zone grids that are connected in a one-to-one, 
patched, or overset manner, and can employ grid sequencing, multigrid, and local time stepping when accelerating 
convergence to steady state. Upwind-biased spatial differencing is used for the inviscid terms, and flux limiting is used 
to obtain smooth solutions in the vicinity of shock waves, when present. Viscous terms are centrally differenced, and 
cross-diffusion terms are neglected. 

The CFL3D code is advanced in time with an implicit approximate factorization method. The implicit derivatives 
are written as spatially first-order accurate, which results in block tridiagonal inversions for each sweep. However, 
for solutions that utilize Roe flux-difference splitting, 19 the block tridiagonal inversions are further simplified using a 
diagonal algorithm with a spectral radius scaling of the viscous terms. 

In time-accurate mode, CFL3D uses pseudo-time stepping with multigrid and achieves second order temporal 
accuracy. With pseudo-time stepping, subiterations are used to reduce the linearization and factorization errors, and 
advance the solution in pseudo-time to the next physical time. For a non-deforming mesh, the time dependent com- 
pressible Navier-Stokes equations can be written as 


1 dQ 
J dt 


= R{Q) 


( 1 ) 


where Q is the vector of conserved variables [p, pu, pv, pw,e] T , J is the Jacobian of the generalized coordinate 
transformation, and R is the right-hand-side: 
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made up of gradients of the inviscid and viscous flux terms. Eq. (1 ) can be discretized with backward differencing: 
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The superscripts n, n+ 1, and n - 1 indicate the time level. When </> = 0 the method is first-order temporally accurate, 
and when </> = 1/2 the method is second-order temporally accurate. To employ pseudo-time stepping, a pseudo time 
term is added to Eq. (3), and the equation is iterated in m, where m is the subiteration counter. After linearizing R and 
performing some additional manipulation, one obtains: 
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where r is the pseudo time parameter, AQ m = Q m+ 1 — Q m , and 

0 = ®^. (5, 

aQ <9Q dQ 

Eq. (4) is approximately factored and written in primitive variable form; it is solved as a series of sweeps in each 
coordinate direction. Additional details are given in Krist et al. 15 

In CFL3D, turbulence models are implemented uncoupled from the mean-flow equations. They are solved using a 
three-factor implicit approximate factorization approach. The advection terms are discretized with first-order upwind 
differencing. The production source term is solved explicitly, while the advection, destruction, and diffusion terms are 
treated implicitly. 

For time-accurate URANS computations, it can be shown 20 that the decomposition for forced periodic cases can 
be written as 

/(x,t) = </(x,t)) +/'(x,t) , (6) 

where /(x, t) is the total signal, the ( ) term is the phase average, and /'(x, t) is the turbulence. This form resembles 
standard Reynolds splitting, except that the flow is split into a slowly varying mean (phase average) and a random 
fluctuating part. As a result, the final conservation equations in terms of phase-averaged variables ( ) are identical in 
form to the standard Reynolds-averaged equations. In other words, when URANS is used for a forced periodic flow 
field such as that produced by a synthetic jet, the computation (depending on the turbulence model and the case) often 
eventually settles down to a nearly exactly repeatable periodic variation. If and when this repeatability is attained, any 
single point during the cycle corresponds to a phase-averaged result from the experiment. 

A schematic diagram showing the hump and the boundary conditions applied is shown in Fig. 1 . The flow control 
slot was located near x/c — 0.65 on the back side of the hump, near where the flow separated in the uncontrolled 
state. A close-up of the region near the slot is shown in Fig. 2. The chord of the hump model was c — 420 mm. The 
freestream Mach number was M — 0.1 (U re f = 34.6 m/s), and the Reynolds number was Re = 936,000 per hump 
chordlength. For the oscillatory case, the oscillation frequency was 138.5 Hz, and the peak velocity out of the slot was 
approximately 26 - 27 m/s. 

All computations performed for this case were 2-D. The grid used was a 2-D structured grid with upper wall 
contouring to approximately account for tunnel blockage. This grid contained 4-zones connected in a 1 -to- 1 fashion, 
and approximately 210,000 total grid points. The internal slot and cavity were included in the grid. Also, a medium- 
level grid was made from the fine grid by extracting every-other point in each coordinate direction. 

Two time steps were investigated: one with 360 time steps per period of the forcing frequency and one with 180 
time steps per period. A total of 20 subiterations were employed per time step. For this case, this number of subitera- 
tions was enough to reduce the Z/ 2 -norm of the subiteration density residual by more than 3 orders of magnitude. See 
Fig. 3, which shows subiteration residual (as well as lift coefficient) history over the course of approximately 7 time 
steps (with 20 subiterations per time step) during part of the unsteady cycle. 

The boundary conditions were as follows. At the floor and hump surfaces, as well as at the side walls inside the 
cavity, solid wall adiabatic boundary conditions were applied. At the front of the grid, which extended to x/c = —6.39, 
far-field Riemann-type boundary conditions were applied. This boundary condition is essentially a non-reflective 
freestream condition; the location x/c = -6.39 was chosen so that the naturally-developing fully-turbulent boundary 
layer in the computation reached the same thickness as the experimentally measured value of approximately 5/c = 
0.073 at the location x/c = —2.14. At the downstream boundary (at x/c = 4.0) the pressure was set at p/p re f = 
0.99962, and all other quantities were extrapolated from the interior of the domain. This back pressure was chosen 
to achieve approximately the correct inflow conditions for steady flow. The top tunnel wall was treated as an inviscid 
wall for all of the computations shown here. Early investigations also explored the effect of modeling the top wall 
as a viscous surface, but this choice made only very small differences in the results. At the bottom of the cavity, the 
boundary condition for the oscillatory case set the velocity components as follows: 

[7 = 0 V = [(/>V% ax / p\cos(2nFt) (7) 

where F is the frequency and t is the time, and ( pV) max was chosen in order to achieve a maximum velocity magnitude 
near to the target of approximately 26 - 27 m/s out of the slot during the cycle. Figure 4 shows the velocity components 
along the slot exit at four different times during the cycle. Maximum blowing occurred at phase = 170°, and maximum 
suction at phase = 350°. The value of (pU) max used to achieve this condition was {pV) max = 0.001p re fa re f, where 
a re f is the reference speed of sound. Along with the above boundary condition given by Eq. (7), the density and 
pressure at the bottom of the cavity were extrapolated from the interior of the domain. Figure 4 also shows the effect 
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of time step on the solution at the slot exit location. There were differences, but these were very small. Although not 
shown, other aspect of the solutions elsewhere in the field using the two different time steps also exhibited only very 
minor differences. This insensitivity to time step at these levels is consistent with results from a different synthetic jet 
study. 20 For all the remaining unsteady results presented herein, the smaller time step corresponding to 360 steps per 
cycle was employed. 

As mentioned earlier, it was discovered during the CFDVAL2004 workshop 2 that the side plates used in the ex- 
periment to improve the flow quality caused blockage that, if not accounted for in the CFD, resulted in computed Cp 
levels that were noticeably different over the hump. This difference was later verified experimentally by obtaining Cp 
data for the no-flow-control case without side plates. One way to approximately account for this blockage in a 2-D 
computation is to modify the upper wall shape to reflect the tunnel area change due to the presence of the plates. This 
was done for the current computations. The differences with and without modified top wall are shown in Fig. 5 for the 
(steady-state) no-flow-control and steady suction cases. In the suction case, note that the large spike in Cp at the slot 
is due to flow acceleration around the upstream edge of the slot lip. 

The effect of grid was investigated for the oscillatory case, and sample results are shown in Fig. 6. Almost no 
differences were visible between long-time-average surface pressures, phase-averaged surface pressures, or velocity 
profiles using the two grids. There were only small differences (generally less than 5% maximum difference) in 
computed turbulent shear stress profiles. Based on this grid study, the medium grid (approximately 53,000 grid points) 
is believed to be fine enough to adequately resolve the numerics of the hump model flow using URANS. Nevertheless, 
for the oscillatory condition the fine grid (approximately 210,000 grid points) was employed to insure that these 
computations possessed more than sufficient grid resolution. 

III. Results 

Final long-time-average mean Cp results for the oscillatory case are shown in Fig. 7. The pressures were predicted 
reasonably well prior to separation, although the computed levels were still as much as 3 - 6% in error over the center 
portion of the hump, despite the use of the modified top wall to account for blockage. Some of this error is likely due 
to the fact that the oscillatory case used a different experimental test configuration underneath the splitter plate, which 
had the overall effect of increasing the suction over the hump model slightly compared to the original configuration. 3 
The use of a lower back pressure in the CFD would improve the correlation. 

The pressure levels in the separated region were missed by a substantially larger amount, and indicated too large 
of a (mean) separation bubble compared with experiment. There were small differences between results from the 
three turbulence models. Table 1 shows measured and predicted separation and reattachment locations (note that only 
the experimental reattachment location for the no-flow-control and steady suction cases were measured accurately 
with oil flow; other experimental values were approximately obtained from particle image velocimetry (PIV) data). 
As shown in the table, CFD predicted separation location fairly accurately (less than 2% error), but overpredicted the 
reattachment location by between 10 - 25%. In fact, the computations predicted only a relatively small change between 
the reattachment location for the no-flow-control condition and the oscillatory condition, whereas the experiment 
indicated that oscillatory control had a fairly significant impact on this aspect of the flowfield. 

The unsteady phase-averaged pressures at four different phases are given in Fig. 8. CFD did an excellent job 
predicting the general character of the unsteady surface pressure coefficients during the forced cycle, including the 
relative effects due to strength and location of the shed vortices that convected downstream. However, the actual Cp 
levels tended to be somewhat low in magnitude (by between 0. 1 - 0.2) for 0.65 < x/c < 0.8, similar to the differences 
seen in the time-averaged results in Fig. 7. The SA model results tended to smear out the vortex more than SST or 
EASM-ko after it passed beyond x/c = 1, in better agreement with the experimental data. 

Phase-averaged u-velocity profiles and turbulent shear stress profiles at three different stations are shown in Figs. 9 
-14. At x/c = 0.66, which was just downstream of the slot, the flow was still attached in the mean (with small 
separation during part of the cycle). During the part of the cycle near phase = 170°, the velocity profiles showed the 
influence of the fluid expulsion from the slot in the form of an inflection in the profile. Velocities near the knee were 
underpredicted at phases 80° and 350°. EASM-ko did the best job overall capturing the peak magnitudes of turbulent 
shear stress at this location. 

At x/c = 0.8, which was in the middle of the separation region, the computations predicted the velocity profiles 
very well in general, including their variation with phase, although there were some specific differences between the 
turbulence models. For example, SST and EASM-ko tended to produce more rounded profiles near the wall than SA. 
All three turbulence models dramatically underpredicted the turbulent shear stress magnitudes at this location within 
the separated region. Probably as a result of this poor prediction of the turbulence, the computed flow field possessed 
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too little mixing and it reattached much later than the experiment. This late reattachment can be seen in Figs. 13 and 
14 at x/c — 1.0. At this location, the experiment has just reattached (in the mean), but the computations remained 
separated over all phases of the cycle. 

Plots showing phase-averaged u - velocity contours at phase = 170° (maximum blowing part of the cycle) are 
given in Fig. 15. At this phase, the experiment exhibited a large area of negative u - velocity over the entire back of 
the hump, with two concentrated regions of higher- velocity fluid: one centered near x/c = 0.72 and a larger one 
centered near x/c = 0.9. These negative velocity regions were a result of convecting vortices (positioned above 
the regions) with clockwise rotation. The computations exhibited qualitatively similar streamwise positioning of the 
reverse-flow concentrated regions, but the negative it-velocity extended further downstream, and the concentrated 
region at x/c = 0.9 was predicted to be larger than the one in the experiment. The SST and EASM-ko results 
exhibited a third concentrated reverse-flow region further downstream near x/c = 1.16 (shed from an earlier cycle), 
whereas SA smeared this feature out. Contours of turbulent shear stress at the same phase are shown in Fig. 16. The 
magnitude of the computed levels in the separated region were too low compared to experiment by a factor of at least 2 
- 3. This figure also shows how SST and EASM-ko tended to maintain more compact regions of turbulence, whereas 
SA diffused the regions more. The long-time-average streamlines are shown in Fig. 17. All three turbulence models 
yielded similar results, predicting larger mean separated regions than experiment. 

Details near the slot are shown for four different phases in Figs. 18-21. The computed velocities were qualitatively 
similar to the experimental results. The separation point oscillated between the slot corner and x/c = 0.665 - 0.666. 
The PIV data was not reliable near the wall, so the precise location of separation in the experiment could not be 
ascertained. However, it appears that CFD may have yielded separation somewhat upstream of experiment in general. 
The three turbulence models predicted similar results for separation location, except at phase = 350° where SST 
and EASM-ko predicted separation slightly downstream of SA. In contrast to the velocity contours, the contours 
of -u'v f indicated dramatic differences between the peak levels given by CFD and experiment. This is consistent 
with the results shown in Fig. 16, and indicates that the turbulence models did not capture the physics correctly even 
immediately downstream of separation. 

In order to determine whether low computed turbulence levels in the separated region can cause delayed reattach- 
ment, a crudely modified SA model was applied to the no-flow-control case. The turbulence model was (arbitrarily) 
modified to produce double the eddy viscosity in an area near the separation bubble, but away from the wall. Resulting 
streamlines from the converged steady state solution are shown in Fig. 22, compared to the standard SA result. (In the 
figure, dashed lines indicate the approximate region in which the eddy viscosity was doubled in the modified case.) 
This ad hoc procedure resulted in a significantly smaller separation bubble, lending support to the conjecture that 
deficiencies in the standard turbulence models are likely to blame for the overpredicted bubble lengths for the hump 
model cases in general. If the turbulence models produced higher magnitude of turbulent shear stress levels in the 
separated region (in better agreement with experiment), then predicted separation extents would decrease, as indicated 
by this simple numerical test. 


IV. Conclusions 

This oscillatory hump case represents a high-quality, nominally two-dimensional experimental dataset that appears 
to expose an inherent flaw in current state-of-the-art RANS and URANS methodologies: the magnitude of the turbulent 
shear stress in the separated region of this configuration is underpredicted, and consequently the separation length is 
overpredicted. As was seen in the CFDVAL2004 workshop, this same flaw was also evident for the hump model in 
steady flow (no-flow-control and steady suction) using a wide variety of codes and turbulence models. The question 
now is whether current turbulence models can be corrected to account for the missed physics, or whether it will be 
necessary to use more computationally expensive, sophisticated methods such as 3-D LES or blended RANS-LES to 
predict these types of synthetic jet flowfields more accurately. 

Overall, 2-D time-accurate URANS showed itself capable of predicting some of the general mean-flow charac- 
ter of this complex unsteady flow, including the relative strength and location of convected vortical flow structures 
as manifested in their effects on time-dependent surface pressures. Three different turbulence models made only 
relatively minor differences in the overall results. However, because of its limitations regarding underprediction of 
turbulent shear stress magnitudes and overprediction of the separation bubble extent, the adequacy of state-of-the-art 
URANS as a design tool for synthetic jets involving separated flow is questionable. Its suitability will likely depend 
on particulars of the flowfield in question, as well as on the level of accuracy required. 
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Table 1. Separation and reattachment locations, x/c 


condition 

model 

sep. loc. 
experiment 

sep. loc. 
CFD 

error 

reatt. loc. 
experiment 

reatt. loc. 
CFD 

error 

no-flow-control 

SA 

« 0.67 

0.661 

1.3% 

1.11 ±0.003 

1.24 

11.7% 

steady suction 

SA 

« 0.68 

0.669 

1.6% 

0.94 ± 0.005 

1.10 

17.0% 

oscillatory (mean) 

SA 

« 0.67 

0.663 

1.0% 

« 0.99 

1.22 

23.2% 

oscillatory (mean) 

SST 

« 0.67 

0.665 

0.7% 

« 0.99 

1.19 

20.2% 

oscillatory (mean) 

EASM-ko 

« 0.67 

0.665 

0.7% 

« 0.99 

1.19 

20.2% 
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Figure 1. Diagram of hump model showing boundary conditions employed. 




x/c x/c 


Figure 2. Close-up of slot exit region; medium-level grid (every other grid point in each coordinate direction) shown. 
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Figure 3. Convergence of subiteration residual during time-accurate oscillatory computation, fine grid, SA. 


thick lines = 360 steps/cycle 
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Figure 4. Velocity components at the slot exit for oscillatory case at 4 different phases in the cycle, fine grid, SA. 
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(a) no-flow-control 




Figure 5. Surface pressure coefficients for no-flow-control and steady suction cases, showing effect of side plate blockage in the experiment 
(no-flow-control case only) and effect of using upper wall contouring to account for blockage in 2-D CFD computations, fine grid, SA. 
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Figure 6. Effect of grid on results for oscillatory case, SA. 
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Figure 7. Long-time-average surface pressure coefficients for oscillatory case. 
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Figure 8. Phase-averaged surface pressure coefficients behind the hump for oscillatory case. 
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Figure 9. Phase-averaged u- velocity profiles at x/c = 0.66 for oscillatory case. 
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Figure 10. Phase-averaged —mV profiles at x/c = 0.66 for oscillatory case. 
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Figure 11. Phase-averaged w-velocity profiles at x/c = 0.8 for oscillatory case. 
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Figure 12. Phase-averaged —u'v' profiles at x/c = 0.8 for oscillatory case. 
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Figure 13. Phase-averaged u-velocity profiles at x/c= 1.0 for oscillatory case. 
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Figure 14. Phase-averaged —u'v' profiles at x/c = 1.0 for oscillatory case. 
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Figure 15. Phase-averaged ir-velocity contours at phase = 170° for oscillatory case; dashed lines indicate regions of zero or negative 
u-velocity (experimental contours under the hump surface should be ignored). 
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Figure 16. Phase-averaged — u'v' contours at phase = 170° for oscillatory case (experimental contours under the hump surface should be 
ignored). 
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Figure 17. Long-time-average streamlines for oscillatory case. 
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Figure 18. Phase-averaged u-velocity and - u'v ' contours near the slot at phase = 80° (separation location for the CFD results are 
indicated by arrows, and white line located y/c = 0.0012 above surface indicates region below which experimental PIV data is considered 
potentially unreliable due to wall reflection). 
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Figure 19. Phase-averaged ir-velocity and —u'v' contours near the slot at phase = 170°, the peak blowing part of the cycle (separation 
location for the CFD results are indicated by arrows, and white line located y/c = 0.0012 above surface indicates region below which 
experimental PIV data is considered potentially unreliable due to wall reflection). 
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Figure 20. Phase-averaged u-velocity and —u'v f contours near the slot at phase = 260° (separation location for the CFD results are 
indicated by arrows, and white line located y/c = 0.0012 above surface indicates region below which experimental PIV data is considered 
potentially unreliable due to wall reflection). 
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Figure 21. Phase-averaged u-velocity and —u'v' contours near the slot at phase = 350°, the peak suction part of the cycle (separation 
location for the CFD results are indicated by arrows, and white line located y/c = 0.0012 above surface indicates region below which 
experimental PIV data is considered potentially unreliable due to wall reflection). 
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Figure 22. Numerical experiment showing streamlines for no-flow-control case, comparing results from standard SA model and SA with 
eddy viscosity arbitrarily doubled in region indicated by dashed lines (starting approximately 0.006 above the wall); this test demonstrates 
how increasing the modeled eddy viscosity in the separated region would decrease the computed bubble size. 
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